{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting run with 1 MPI rank(s) at : 2020-08-18 10:41:49.147374\n"
]
}
],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"\n",
"from triqs.plot.mpl_interface import plt\n",
"from triqs_tprf.tight_binding import create_square_lattice\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"plt.style.use('notebook.mplstyle')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linearized Eliashberg equation on the attractive Hubbard model\n",
"\n",
"A simple example for the usage of the linearized Eliashberg equation is the attractive Hubbard model on a square lattice.\n",
"It is not only fast and simple to setup, but the particle-hole symmetry of the Hubbard model also serves as a benchmark for the correctness of such a calculation.\n",
"\n",
"In the following we will first introduce the Hubbard model and present the semi particle-hole transformation and its usage.\n",
"Then we will show how to solve the linearized Eliashberg equation for this model to find the superconducting phase transition using TRIQS and TPRF routines.\n",
"\n",
"If you want a more detailed study of this problem, checkout this [notebook](https://github.com/TRIQS/tprf/blob/eliashberg/benchmark/eliashberg/particle_hole_transformation/PHT_Hubbard_Model.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hubbard model on a square lattice\n",
"\n",
"The particle-hole symmetric Hamiltonian for the Hubbard model is given by\n",
"\n",
"$$\n",
"H=-t \\sum_{\\langle j, l\\rangle \\sigma}\\left(c_{j \\sigma}^{\\dagger} c_{l \\sigma}+c_{l \\sigma}^{\\dagger} c_{j \\sigma}\\right)+U \\sum_{j}\\left(n_{j \\uparrow}-\\frac{1}{2}\\right)\\left(n_{j \\downarrow}-\\frac{1}{2}\\right)-\\mu \\sum_{j}\\left(n_{j \\uparrow}+n_{j \\downarrow}\\right)\\,,\n",
"$$\n",
"\n",
"here $c_{j\\sigma}^{\\dagger}$ creates an electron on site $j$ with spin $\\sigma$ while $c_{j\\sigma}$ destroys such an electron, further the operator $n_{j\\sigma}$ count the number of electrons on site $j$ with spin $\\sigma$.\n",
"\n",
"The first term describes the kinetic energy of the electrons, which can be interpreted as an electron with spin $\\sigma$ *hopping* from site $l$ to site $j$ and vice versa.\n",
"Here the angular braket under the sum means that we only take *hopping* terms between neighboring lattice sites into account and the energy that is gained by such a *hopping* process is given by $t$.\n",
"\n",
"The second term describes the repulsive interaction between the electrons.\n",
"This repulsion is crudely approximated in the Hubbard model in the sense, that electrons only *see* each other if they occupy the same lattice site.\n",
"The energy that is needed to have a lattice site doubly occupied is given by $U$.\n",
"\n",
"The last term describes the filling of the lattice via an energy offset by the chemical potential $\\mu$.\n",
"\n",
"The Hubbard model is therefore defined by the parameters $t$, $U$ and $\\mu$, but we also need to know the temperature $T$ at which we shall observe the Hubbard model.\n",
"We will record these parameters using the `ParameterCollection` class of `triqs_tprf.ParameterCollection` module."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"T = 1000\n",
"U = 1.0\n",
"mu = 0.0\n",
"nk = 32\n",
"norb = 1\n",
"nw = 50\n",
"spin = False\n",
"t = 1.0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from triqs_tprf.ParameterCollection import ParameterCollection\n",
"\n",
"hubbard = ParameterCollection( # -- Model Parameter\n",
" norb=1, # Number of orbitals.\n",
" t=1.0, # Hopping to nearest neighbor\n",
" U=1.0, # Strength of the on-site interaction\n",
" mu=0.0, # Chemical potential determining the filling.\n",
" T=1000, # Temperature.\n",
" spin=False, # Treat indices only for orbital character.\n",
" \n",
" # -- Technical parameter\n",
" nk=32, # Number of points in one dimension considered in the Brillouin zone.\n",
" nw=50, # Number of Matsubara points in positive dimension.\n",
" )\n",
"hubbard"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Semi particle-hole transformation\n",
"\n",
"The Hubbard model with only nearest neighbor hopping inhibts a useful semi particle-hole symmetry, which we will present without any rigorous derivations.\n",
"\n",
"First off, a bipartite lattice can be subdivided into two sublattices for which every lattice site on one of them only has neighboring sites from the other sublattice. This is the case for our square lattice.\n",
"With this knowledge we can define the particle-hole transformation (PHT) for the creation and annihilation operators as \n",
"\n",
"$$\n",
"c^\\dagger_{j \\sigma} \\xrightarrow{\\mathrm{PHT}} (-1)^{j} d_{j \\sigma} \\,,\\quad\n",
"\\mathrm{and}\\quad\n",
"c_{j \\sigma} \\xrightarrow{\\mathrm{PHT}}(-1)^{j}d_{j \\sigma}^{\\dagger}\\,,\\\\\n",
"$$\n",
"\n",
"where $j$ is either $0$ for one sublattice and $1$ for the other.\n",
"\n",
"If we use the PHT only on one spin specices it is called a semi particle-hole transformation (SPHT).\n",
"This means\n",
"\n",
"$$\n",
"c_{j \\uparrow}^{\\dagger} \\xrightarrow{\\mathrm{SPHT}} d_{j \\uparrow}^{\\dagger}\\quad\\mathrm{and}\\quad\n",
"c_{j \\uparrow} \\xrightarrow{\\mathrm{SPHT}} d_{j \\uparrow}\\\\\n",
"c_{j \\downarrow}^{\\dagger} \\xrightarrow{\\mathrm{SPHT}} (-1)^j d_{j \\downarrow}\\quad\\mathrm{and}\\quad\n",
"c_{j \\downarrow} \\xrightarrow{\\mathrm{SPHT}} (-1)^j d_{j \\downarrow}^{\\dagger}\\,,\n",
"$$\n",
"\n",
"and for the number operators\n",
"\n",
"$$\n",
"n_{j \\uparrow}\\xrightarrow{\\mathrm{SPHT}} \\tilde{n}_{j \\uparrow}\n",
"\\quad\\mathrm{and}\\quad\n",
"n_{j \\downarrow} \\xrightarrow{\\mathrm{SPHT}} 1-\\tilde{n}_{j \\downarrow}\\,.\n",
"$$\n",
"\n",
"One can convince themselves, that the kinetic part of the Hubbard model on a square lattice with only nearest neighbor hopping is invariant under a SPHT.\n",
"The interaction term on the other hand is not and changes sign\n",
"\n",
"$$\n",
"\\left(n_{j \\uparrow}-\\frac{1}{2}\\right)\\left(n_{j \\downarrow}-\\frac{1}{2}\\right) \\xrightarrow{\\mathrm{SPHT}}\n",
"-\\left(\\tilde{n}_{j \\uparrow} - \\frac{1}{2}\\right)\\left(\\tilde{n}_{j \\downarrow}-\\frac{1}{2}\\right)\\,.\n",
"$$\n",
"\n",
"Therefore the SPHT maps the repulsive Hubbard model with $U$ to the attractive one with $-U$.\n",
"Additionally to that, if our Hubbard model is not a half-filling, the chemical potential term is also not invariant under a SPHT\n",
"\n",
"$$\n",
"n_{j \\uparrow}+n_{j \\downarrow} \\xrightarrow{\\mathrm{SPHT}}\n",
"1 + \\left(\\tilde{n}_{j \\uparrow}-\\tilde{n}_{j \\downarrow}\\right)\\,,\n",
"$$\n",
"\n",
"and transforms into a Zeeman term.\n",
"\n",
"To summarize the SPHT maps the Hubbard Hamiltonian with interaction strength $U$ and chemical potential $\\mu$ to a Hubbard Hamiltonian with interaction strength $-U$, a chemical potential of $0$ and an additional Zeeman term of strength $\\mu$.\n",
"If we therefore calculate an observable $A$ in the attractive Hubbard model we know the observable $B\\xleftarrow{\\mathrm{SPHT}}A$ in the repulsive one.\n",
"\n",
"For example, the in-plane antiferromagnetic (AFM) phase of the repulsive model is connected to the superconducting (SC) phase in the attractive one, as can be seen by using the SPHT on the ladder operators of the spin\n",
"\n",
"$$\n",
"S^+_j = S^x_j + iS^y_j = c_{j \\uparrow}^{\\dagger}c_{j \\downarrow} \\xrightarrow{\\mathrm{SPHT}} (-1)^j d_{j \\uparrow}^{\\dagger}d_{j \\downarrow}^{\\dagger} = \\tilde{\\Delta}^{\\dagger}\\,,\\\\\n",
"S^-_j = S^x_j - iS^y_j = c_{j \\downarrow}^{\\dagger}c_{j \\uparrow} \\xrightarrow{\\mathrm{SPHT}} (-1)^j d_{j \\downarrow}d_{j \\uparrow} = \\tilde{\\Delta}\\,.\n",
"$$\n",
"\n",
"The $x$- and $y$-components of the spin operator are transformed to the complex superconducting oder parameter with a phase factor.\n",
"This means, that if we find a staggered in-plane spin phase in the repulsive model at some $U$, we will see a homogeneous superconducting phase at $-U$.\n",
"\n",
"The following plot shows this symmetry in a T-U phase diagram, which was previously calculated [here](https://github.com/TRIQS/tprf/blob/eliashberg/benchmark/eliashberg/particle_hole_transformation/PHT_Hubbard_Model.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import SVG\n",
"\n",
"SVG(filename='./plots/SPHT_hubbard_phase_diagram.svg')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The right hand side of the phase diagram was calculated using an attractive Hubbard model at half-filling, i.e. $\\mu=0.0$, with an additional Zeeman term of strength $\\xi$.\n",
"The phase boundary was then calculated using the random phase approximation (RPA), i.e. using a frequency independent and local vertex $U$, to calculate $\\langle S^xS^x \\rangle$ and increasing $U$ until divergence.\n",
"For details on how to use TPRF for that see this [tutorial](Square lattice susceptibility.ipynb#Random-phase-approximation-(RPA)).\n",
"\n",
"The left hand side of the phase diagram was calculated using an repulsive Hubbard model with $\\mu=\\xi$, i.e. the SPHT mapping of the right hand side model.\n",
"Then the linearized Eliashberg equation was solved for this model for various interaction strengths until the superconducting phase transition was found."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving the linearized Eliashberg equation\n",
"\n",
"The linearized Eliashberg equation is given in a very simplified form by\n",
"\n",
"$$\n",
"\\lambda\\Delta = \\Lambda \\Delta\\,,\n",
"$$\n",
"\n",
"where $\\Delta$ is a gap function and $\\lambda$ the largest eigenvalue of the matrix $\\Lambda$.\n",
"Here $\\Lambda$ is the product of the particle-particle vertex $\\Gamma$ and the Green's function $G$.\n",
"If $\\lambda=1$ a phase transition to a superconducting state is found.\n",
"For further information see the documentation [here](../theory/eliashberg.rst).\n",
"\n",
"To solve it in the same order as the RPA we need the non-interacting Green's function $G^{(0)}$ and the (singlet) particle-particle vertex $\\Gamma$ is approximated by a fequency independent and local vertex $2U$.\n",
"\n",
"First we setup our model by using previously established `ParameterCollection` `hubbard` as a template and alter it to the parameters of the repulsive Hubbard model.\n",
"To do this use its method `alter` and supply the parameters that shall be changed as keywords.\n",
"We set the chemical potential to $\\mu=\\xi$ and the interaction strength to $U=-1.41$, which is close to the superconducting phase boundary for $T=1000\\,\\mathrm{K}$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"T = 1000\n",
"U = -1.41\n",
"mu = 0.1\n",
"nk = 32\n",
"norb = 1\n",
"nw = 50\n",
"spin = False\n",
"t = 1.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xi = 0.1\n",
"\n",
"repl_hubbard = hubbard.alter(mu=xi, U=-1.41)\n",
"repl_hubbard"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To construct a representation of the kinetic part of the Hubbard model use the `create_square_lattice` function of the `triqs_tprf.tight_binding` module.\n",
"Give it as keywords the number of orbitals `norb` and the hopping energy `t`, which can comfortably be accessed from `repl_hubbard`.\n",
"It returns a `TBLattice` object."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"triqs_tprf.tight_binding.TBLattice"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from triqs_tprf.tight_binding import create_square_lattice\n",
"\n",
"H = create_square_lattice(norb=repl_hubbard.norb, t=repl_hubbard.t)\n",
"type(H)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To get the dispersion relation, use its member function `fourier` and enter a mesh obtained using the other member function `get_kmesh` on the Brillouin zone as a tuple.\n",
"It returns the dispersion relation stored as a `Gf` object."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"triqs.gf.gf.Gf"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e_k = H.fourier(H.get_kmesh(n_k=(repl_hubbard.nk, repl_hubbard.nk, 1)))\n",
"type(e_k)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the dispersion relation we can construct the non-interacting Green's function.\n",
"To do this first create a `MeshImFreq` object with a fermionic statistic and the wished inverse temperature $\\beta$ and number of points.\n",
"Then use the `lattice_dyson_g0_wk` function of the `triqs_tprf.lattice` module and supply it with the dispersion relation, the `MeshImFreq` object and a chemical potential.\n",
"This yields the non-interacting Green's function as a `Gf` object.\n",
"(You can use the `temperature_to_beta` function of the `triqs_tprf.utilities` to calculate the inverse temperature in $\\mathrm{eV}$ from the temperature in $\\mathrm{Kelvin}$.)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"triqs.gf.gf.Gf"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from triqs_tprf.lattice import lattice_dyson_g0_wk\n",
"from triqs.gf import MeshImFreq\n",
"from triqs_tprf.utilities import temperature_to_beta\n",
"\n",
"beta = temperature_to_beta(repl_hubbard.T)\n",
"wmesh = MeshImFreq(beta=beta, S='Fermion', n_max=repl_hubbard.nw)\n",
"g0_wk = lattice_dyson_g0_wk(mu=repl_hubbard.mu, e_k=e_k, mesh=wmesh)\n",
"type(g0_wk)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The particle-particle vertex $\\Gamma$ in this case must be manually created.\n",
"To do this first create a `MeshImFreq` in the same manner as before, but this time with a bosonic statistic.\n",
"Then use this `MeshImFreq` and combine it with the momentum mesh of the non-interacting Green's function to create a `MeshProduct` object.\n",
"Use this `MeshProduct` object to create a `Gf` object with a `target_shape` that corresponds to a two-particle object, e.g. `(1, 1, 1, 1)`.\n",
"Then overwrite its data to be constant to $U$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"from triqs.gf import Gf, MeshProduct\n",
"\n",
"wmesh_boson = MeshImFreq(beta=temperature_to_beta(repl_hubbard.T), S='Boson', n_max=repl_hubbard.nw)\n",
"wmesh_boson_kmesh = MeshProduct(wmesh_boson, g0_wk.mesh[1])\n",
"gamma_pp = Gf(mesh=wmesh_boson_kmesh, target_shape=g0_wk.target_shape*2)\n",
"gamma_pp.data[:] = 2*repl_hubbard.U"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With this we have all ingredients for the linearized Eliashberg equation.\n",
"To solve it use the `solve_eliashberg` function of the `triqs_tprf.eliashberg` module and supply it with the particle-particle vertex and the non-interacting Green's function.\n",
"It returns the eigenvalues and eigenvectors, i.e. the gap functions.\n",
"\n",
"The returned maximum eigenvalue is $\\lambda \\approx 1$ as expected, because we are close to the phase boundary."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.9999705817306739"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from triqs_tprf.eliashberg import solve_eliashberg\n",
"\n",
"Es, eigen_modes = solve_eliashberg(gamma_pp, g0_wk)\n",
"Es[0]"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}